const FastMath = require('./FastMath');

function LUDecomposition(matrix, singularityThreshold) {
  var DEFAULT_TOO_SMALL = 1e-11;

  singularityThreshold = singularityThreshold || DEFAULT_TOO_SMALL;

  if (!matrix.isSquare()) {
    throw new NonSquareMatrixException(matrix.getRowDimension(),
      matrix.getColumnDimension());
  }

  var m = matrix.getColumnDimension();
  this.lu = matrix.getData();
  this.pivot = new Array(m);
  this.cachedL = null;
  this.cachedU = null;
  this.cachedP = null;

  // Initialize permutation array and parity
  for (var row = 0; row < m; row++) {
    this.pivot[row] = row;
  }
  this.even = true;
  this.singular = false;

  // Loop over columns
  for (var col = 0; col < m; col++) {

    // upper
    for (var row = 0; row < col; row++) {
      var luRow = this.lu[row];
      var sum = luRow[col];
      for (var i = 0; i < row; i++) {
        sum -= luRow[i] * this.lu[i][col];
      }
      luRow[col] = sum;
    }

    // lower
    var max = col; // permutation row
    var largest = -Infinity;
    for (var row = col; row < m; row++) {
      var luRow = this.lu[row];
      var sum = luRow[col];
      for (var i = 0; i < col; i++) {
        sum -= luRow[i] * this.lu[i][col];
      }
      luRow[col] = sum;

      // maintain best permutation choice
      if (FastMath.abs(sum) > largest) {
        largest = FastMath.abs(sum);
        max = row;
      }
    }

    // Singularity check
    if (FastMath.abs(this.lu[max][col]) < singularityThreshold) {
      singular = true;
      return;
    }

    // Pivot if necessary
    if (max != col) {
      var tmp = 0;
      var luMax = this.lu[max];
      var luCol = this.lu[col];
      for (var i = 0; i < m; i++) {
        tmp = luMax[i];
        luMax[i] = luCol[i];
        luCol[i] = tmp;
      }
      var temp = this.pivot[max];
      this.pivot[max] = this.pivot[col];
      this.pivot[col] = temp;
      this.even = !this.even;
    }

    // Divide the lower elements by the "winning" diagonal elt.
    var luDiag = this.lu[col][col];
    for (var row = col + 1; row < m; row++) {
      this.lu[row][col] /= luDiag;
    }
  }

  this.getSolver = function () {
    return new Solver(this.lu, this.pivot, this.singular);
  }
}

function Solver(lu, pivot, singular) {
  this.lu = lu;
  this.pivot = pivot;
  this.singular = singular;

  this.getInverse = function() {
    return this.solve(require('./MatrixUtils').createRealIdentityMatrix(pivot.length));
  }

  this.solve = function(b) {
    var m = this.pivot.length;
    if (b.getRowDimension() != m) {
      throw new DimensionMismatchException(b.getRowDimension(), m);
    }
    if (this.singular) {
      throw new SingularMatrixException();
    }
    var nColB = b.getColumnDimension();

    // Apply permutations to b
    var bp = new Array(m);
    for (var row = 0; row < m; row++) {
      bp[row] = new Array(nColB);
      var bpRow = bp[row];
      var pRow = this.pivot[row];
      for (var col = 0; col < nColB; col++) {
        bpRow[col] = b.getEntry(pRow, col);
      }
    }

    // Solve LY = b
    for (var col = 0; col < m; col++) {
      var bpCol = bp[col];
      for (var i = col + 1; i < m; i++) {
        var bpI = bp[i];
        var luICol = lu[i][col];
        for (var j = 0; j < nColB; j++) {
          bpI[j] -= bpCol[j] * luICol;
        }
      }
    }

    // Solve UX = Y
    for (var col = m - 1; col >= 0; col--) {
      var bpCol = bp[col];
      var luDiag = lu[col][col];
      for (var j = 0; j < nColB; j++) {
        bpCol[j] /= luDiag;
      }
      for (var i = 0; i < col; i++) {
        var bpI = bp[i];
        var luICol = lu[i][col];
        for (var j = 0; j < nColB; j++) {
          bpI[j] -= bpCol[j] * luICol;
        }
      }
    }

    const Array2DRowRealMatrix = require('./Array2DRowRealMatrix');
    return new Array2DRowRealMatrix(bp, false);
  }

}

module.exports = LUDecomposition
